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Abstract — Simulating the Space Radiation environment 
with Monte Carlo Codes, such as FLUKA, requires the 
ability to model the interactions of heavy ions as they 
penetrate spacecraft and crew member's bodies. Monte- 
Carlo-type transport codes use total interaction cross 
sections to determine probabilistically when a particular 
type of interaction has occurred. Then, at that point, a 
distinct event generator is employed to determine 
separately the results of that interaction. The space 
radiation environment contains a full spectrum of radiation 
types, including relativistic nuclei, which are the most 
important component for the evaluation of crew doses. 
Interactions between incident protons with target nuclei in 
the spacecraft materials and crew member’s bodies are well 
understood. However, the situation is substantially less 
comfortable for incident heavier nuclei (heavy ions). We 
have been engaged in developing several related heavy ion 
interaction models based on a Quantum Molecular 
Dynamics-type approach for energies up through about 5 
GeV per nucleon (GeV/A) as part of a NASA Consortium 
that includes a parallel program of cross section 
measurements to guide and verify this code development. 
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1. Introduction 

The space radiation environment is complex and 
impossible to re-create in the laboratory. While individual 
components can be provided in ground-level experiments it 
is beyond the current technology to reproduce the combined 
spectra of all the constituents. As such, we are constrained 
to rely on computer simulations to predict the nature of the 
radiation field behind various shielding scenarios that are 
exposed to the space radiation environment. One of the 
computational techniques employed is the so-called Monte- 
Carlo transport method. Based on a representation of the 
geometry and composition of the shielding materials, 
individual incident particles are transported and 
interactions are modeled probabilistically based on the 
known cross sections for them within each particular kind 
of medium. The computer codes themselves generally split 
the calculations into two distinct parts, one of which uses 
the relevant total cross section to ask if any interaction 
occurs, and the second, called an event generator is then 
employed to simulate the details of the interaction 
whenever one takes place. Typical uses call for the 


simulation of very large numbers of incident particles in 
order to reduce the statistical uncertainty in the final 
predictions. Of course, statistics cannot overcome any 
weaknesses in the accuracy of the basic event generators. 

Excellent event generators exist for most kinds of 
interactions, such as those which occur when electrons or 
photons are being propagated, or when individual nuclear 
particles such as protons, neutrons or pions are involved. 
However, the event generators for nucleus-nucleus 
collisions, especially in the energy regime from the 
threshold for such interactions starting around a few 
MeV/A up though about 5 GeV/A are on much less firm 
footing. Part of the reason for this is the wide range of 
physics that occurs in this part of the energy spectrum and 
the involvement of the details of nuclear structure. At 
much higher energies, the situation actually becomes easier 
to model as far as the aspects of interest in predicting dose- 
type effects are concerned, and again adequate event 
generators for those interactions also exist. 

In order to address this need, NASA has formed two related 
Consortia, one devoted to the modeling necessary to 
produce these event generators, and the other charged with 
measuring the cross sections necessary to guide and 
validate the models. [1] One of the reasons that the existing 
models are not as far along as those for other types of 
interactions is indeed the lack of data of the kind that are 
necessary to fine- tune and gain confidence in the models. 
It is not possible to make a detailed enough series of 
measurements with all projectile-target-energy 
combinations to use empirical models that simply fit the 
data. Rather, one must develop physics-based models and 
gain sufficient confidence that they do represent the 
fundamental physics processes correctly enough to be valid 
within reasonable limits when they are called upon to 
simulate interactions for which we have no data. 

To that end, we have selected the RQMD models of Sorge 
as a starting point. [2] One of the most successful 
approaches to modeling the interactions of individual 
particles with nuclei has been the Intra-Nuclear Cascade 
(INC) technique. [3] In this model the target nucleus is 
represented as a collection of protons and neutrons with 
appropriate energy distributions held in a global nuclear 
field. The incident particle is then transported through this 
representation probabilistically using the free-particle 
interaction cross sections with the protons and neutrons, 
and employing a separate model for the interaction of the 
particle with the nuclear field. In these models, a constant 
mean field is assumed throughout. This “mean-field” 
approximation works well when the incident particle is a 
proton and the target is a heavy nucleus. However, when 
the projectile and target are both composite nuclei, this 
mean-field approximation is likely to be less accurate in its 
representation of the situation. 
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Chemists have faced a similar problem in modeling 
molecular dynamics. When a molecule is isolated, its 
structure is typically easily modeled. However, in the 
proximity of other molecules, the atoms in the first 
molecule begin to feel the separate influence of the charges 
of other nearby atoms in adjacent molecules, distorting the 
real field that they see. To overcome this, an iterative 
extension called QMD (Quantum Molecular Dynamics) 
was introduced. [4] That approach actually recalculates the 
field at the end of each small time step felt by every particle 
due to the presence of all nearby atoms, and not just the 
ones in its own molecule. RQMD is an adaptation of that 
technique to the world of nuclear interactions. As in the 
case of the INC models, the nuclei are represented by 
collections of protons and neutrons that are being held 
together by their own mutual nuclear field. The difference 
is that instead of using the mean field method, RQMD 
follows the lead of the chemists and recalculates the field at 
each time step. 

We have a version of the FLUKA transport code that has 
had a modified version of he RQMD code embedded in it to 
handle the nucleus-nucleus collisions up to an energy of 5 
GeV/A.[5] To cover energies above this value we embed 
the DPMJET III event generator. [6] We are also working 
on an extension to an additional event generator to be 
employed at energies below 100 MeV/A based on the 
Boltzmann Master Equation approach. [7] The present 
paper will focus on our progress and some of our results on 
the RQMD event generator and our efforts to come up with 
an improved version as well as an alternate event generator 
using a somewhat modified approach. 


2. The RQMD Event Generator 

The version of RQMD that is included in the present 
release of FLUKA is a modified version of the well known 
code developed by Sorge.[2] This code has an option to run 
in an INC-like mode to reduce computation time in 
instances where the full accuracy is not required. The 
fundamental code has not been changed, but some 
modifications to the outputs in the form of energy re- 
normalizations were required to insure absolute energy 
conservation, which is an essential requirement of the 
overlying FLUKA transport code. 

The RQMD code itself can in principle reach the final state 
directly, but allowing it to do so would be prohibitive in 
execution time. In addition, the approximations used to 
recalculate the nuclear field, which are also necessary to 
keep execution time within reason, do not give a 
particularly impressive result for representing some of the 
correct details of the nuclear structure. As such, the 
strategy is to allow the code to evolve the interaction to 
some intermediate point where essentially all of the hard 


collisions by the constituent particles are over, and then to 
halt the evolution. A separate routine then surveys this 
interim state and attempts to collect the nucleons into 
fragments based on their proximity in phase space. Lastly, 
the remaining particles are placed in the final state event 
generator output buffer, and the trial fragments identified 
in the previous step are run through a pre-equilibrium 
processor and an evaporation code that allows them to de- 
excite into a surviving fragment, which may have boiled-off 
some of its original constituents in the evaporation process. 
The selection of the exact point in the RQMD evolution at 
which the process is turned over to these final state 
handlers is a potentially tunable parameter that can be used 
as data become available to attempt to optimize the fit. 

Figures 1-6 show some results from the RQMD event 
generator as it is exercised within the FLUKA code. These 
fluence plots (number of particles cm' 2 ) were all prepared to 
aid in the planning for proposed future measurements at the 
Alternating Gradient Synchrotron (AGS) at the 
Brookhaven National Laboratory on Long Island. They all 
depict the laboratory scattering angles in annular angle bins 
of one degree for protons above a kinetic energy of 100 
MeV in the laboratory, for all charged pions above 50 MeV 
in the laboratory and for all light ions (Deuterons through 
Boron) above 100 MeV/A. The ordinates are all 
normalized to the yields from a 1 kHz beam incident on a 
one interaction length target. 

All of the figures show a similar behavior with a proton 
distribution peaking near zero (note that the integration 
over the annular area of each angle bin distorts the areal 
density distribution of the tracks, which is centrally 
peaked). The pion scattering angle distribution rises from 
low values at small angles, becoming comparable to or 
exceeding the proton fluences between 15 and 20 degrees. 

One of the important differences in the calculations for 5 
GeV/A incident projectiles with respect to similar 3 GeV/A 
interactions is the relative height of the pion distributions at 
these larger angles. 
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Au ■ Si @ 3.0 GeV/A 



Figure 1 - Laboratory Scattering Angle Distributions for 
protons, pions and light ions (D-B) per second in one 
degree annular bins from a 1 kHz 3.0 GeV/A Gold (Au) 
beam incident on a one interaction length thick Silicon (Si) 
target. 



Figure 2 - This should be compared with Figure 1 . It 
shows the Laboratory Scattering Angle Distributions for 
protons, pions and light ions (D-B) per second in one 
degree annular bins from a 1 kHz 5.0 GeV/A Gold (Au) 
beam incident on a one interaction length thick Silicon (Si) 
target. Note the dominance of the pions starting at about 
13 degrees. 



Figure 3 - Laboratory Scattering Angle Distributions for 
protons, pions and light ions (D-B) per second in one 
degree annular bins from a 1 kHz 3.0 GeV/A Iron (Fe) 
beam incident on a one interaction length thick Iron (Fe) 
target. Note the greater angular range for the light ions 
and the flatter proton distrubution as compared with the Au- 
Si plots in Figures 1 and 2. 

In comparing these plots it is important to consider the 
absolute values of the fluences as a function of both energy 
and the size of the projectile and target. Generally, the 
fluences are bivariant and scale with both of these 
parameters. The plots tend to include mostly the 
contributions from the projectiles only, as the fragments 
from the target typically have laboratory kinetic energies 
below the cutoffs used here. There are differences in the 
plots in terms of the relative ratios of protons to pions and 
in 



Figure 4 - Laboratory Scattering Angle Distributions for 
protons, pions and light ions (D-B) per second in one 
degree annular bins from a 1 kHz 5.0 GeV/A Iron (Fe) 
beam incident on a one interaction length thick Iron (Fe) 
target. Note the similar shape, but at significantly lower 
fluences than in the Au-Si case at 5.0 GeV/A in Figure 2. 
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Figure 5 - Laboratory Scattering Angle Distributions for 
protons, pions and light ions (D-B) per second in one 
degree annular bins from a 1 kHz 3.0 GeV/A Silicon (Si) 
beam incident on a one interaction length thick Silicon (Si) 
target. Compare this with the plot in Figure 3, noting the 
similarity in the shape of the distributions but the difference 
in the absolute magnitudes of the fluences. Also note the 
differences in the light ion fluences in the central most 
angle bin. 

the angular distributions of the light ions, especially at the 
more central angle bins. Also, in Figures 1 and 6, the 
collision is between decidedly asymmetrical projectile and 
target masses. These plots show some of the differences 
one can expect from asymmetrical situations such as given 
here. 



Figure 6 - Laboratory Scattering Angle Distributions for 
protons, pions and light ions (D-B) per second in one 
degree annular bins from a 1 kHz 5.0 GeV/A Carbon (C) 
beam incident on a one interaction length thick Silicon (Fe) 
target. Note that this represents a case of very non- 
symmetric collision participants, and the effects are 
noticeable as compared to the more symmetric participants 
in Figures 2-5. 


3. Improved Event Generators 

We have taken three paths towards improving this existing 
RQMD event generator. The first is to tune up the 
supporting software that does the fragment collecting after 
the event generator evolution is stopped and the 
evaporation of those fragments occurs. The second is to 
delve into the core of the event generator to try and clean 
up the approximations that were made by the original 
author and lead to the energy non-conservation behavior. 
Both of these two approaches keep the same basic code in 
place. The third path is to develop an entirely parallel 


approach, termed the Hamiltonian Molecular Dynamics 
model (HMD) that focuses on the details of the needed 
relativistic transformations internally within the evolution 
processor. None of these efforts is ready to produce results 
for comparison with the existing generator or with data, but 
hopefully by the time of the oral presentation, some 
preliminary results will be forthcoming. We limit ourselves 
here to providing a brief description of the work on this 
third path. 

The HMD model is based on the concept of the constrained 
Hamiltonian as introduced by Dirac, [8] and developed by 
Todorov and Komar. [9, 10] Effective nucleon-nucleon 
potentials are given in a non-relativistic form (e.g. as a 
function of 3-distance). Clearly the challenge is to cast 
nucleon (and thus nuclear) equations of motion which are 
nevertheless manifestly covariant. In moving from a 
Wigner-picture 6-N dimensional quasi-phase space to 
Minkowski space, we see that the number of degrees of 
freedom increase by 2N. Thus, 2N equations of constraint 
allow us to define the world-lines for each particle in the 
system uniquely. 

Potentials are included by making the philosophical 
approximation that the effect of the potentials on system 
trajectory is to alter the transverse momentum distributions, 
so that the argument of the potential functional is 
transformed to the transverse 4-distance, making the 
potentials by definition Poincare scalars. Alternatively, one 
can also show by solution of the system Hamiltonian 
characteristic equations that this choice of argument is the 
single value allowing the Hamiltonian to reach the local 
extremum value required by the variational principle. 

Working according to the constrained Hamiltonian 
formalism, the Hamiltonian is written as a linear 
combination of the (scalar) constraint equations, the 
coefficients being elements from the inverse of the Poisson 
matrix. We believe that the energy conservation issues 
realized by previous authors are overcome via a tedious, but 
necessary individual calculation and combination of single- 
particle Hamiltonians and correlation (potentials and 
Hilbert coherence) energies. This has been shown correct 
for a limited range of systems, and a more general 
validation of this contention is currently ongoing. [2, 11, 12, 
13] 

In its current form, this model makes use of the free- 
particle nucleon-nucleon elastic cross sections to describe 
successive particle collisions, with no inclusion of inelastic 
channel probabilities. We recognize that this is insufficient 
for a full description of the system, both in terms of the 
inclusion of inelastic reactions, and the alteration to some 
form of in-medium cross sections. This addition will be 
made and the current version should be viewed as an 
interim developmental step in the complete process. 
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Because of the reliance of this method on Hamilton- Jacobi 
theory, Pauli blocking, and fragment coalescence, quantum 
final state interactions are included a posteriori rather than 
being a natural consequence of the central theory. Pauli 
blocking is included such that a particle cannot scatter back 
into a state with momentum below the Fermi level in its 
original nucleus. Coalescence is achieved through an again 
tedious, but necessary calculation of combinatorial 
possibilities based on proximity with regard to spatial 
coordinate and momentum, effectively working towards a 
minimization of “coalescence volume”. Final state 
interactions are included by allowing further collisions and 
interactions to occur after the initial set of reactions 
characteristic of the two nuclei passing through the overlap 
or “fireball” region. This includes interactions between 
product pre-fragments. 


4. Conclusions 

We have succeeded in embedding a heavy ion event 
generator based on the RQMD approach into the FLUKA 
Monte-Carlo Transport code. While that event generator 
has some potential shortcomings, it represents one of the 
best alternatives for modeling nuclear interactions that 
occur in collisions of heavy ions at energies from -100 
MeV/A up through -5 GeV/A. This effort was motivated 
principally by NASA’s need to be able to simulate the 
effects of the exposure of spacecraft and crew members to 
the space radiation environment. At present, NASA has 
created a pair of coupled consortia that are charged with 
developing the kind of event generators that have been 
discussed here on the one hand, and with making the 
measurements needed to evaluate the efficacy of those event 
generators on the other. It is NASA’s intention that these 
two consortia function in an iterative matter to increase the 
accuracy of the event generator models to the point where 
the fluences can be predicted to an accuracy of within 25% 
from the simulations. 

The uses that NASA intends to make of these codes include 
the direct aiding in the design of the spacecraft and 
structures that future astronauts will employ in the effort to 
return to the Moon and to explore Mars. These include the 
effects of the Martian atmosphere and the Martian and 
Lunar regoliths as well as for any spacecraft and space suit 
configurations. 

These codes represent the best alternative in the evaluation 
of the radiation field where complex geometries and mixed 
elemental materials are concerned. They nicely 
complement simpler approaches used in estimating this 
same thing and can act as a validation tool for those 
calculations as well as for detailed ground and space-based 
measurements. 


The work on these event generators is far from being 
finished. The greatest need is for accurate data sets over a 
range of projectile, target and energy values so that the 
event generators can be improved and ultimately validated. 
In preparation for these measurements, we are engaged in 
the development of extensions to the present theoretical 
approaches to the physics modeling. We are involved in 
the development of a new HMD event generator that will 
hopefully be ready for testing as the next round of data 
becomes available. 
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